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ABSTRACT 

I utilize the Petrov-Galerkin formulation and develop a new method for solving the unsteady coUi- 
sionless Boltzmann equation in both the linear and nonlinear regimes. In the first order approximation, 
the method reduces to a linear eigenvalue problem which is solved using standard numerical meth- 
ods. I apply the method to the dynamics of a model stellar disk which is embedded in the held of 
a soft-centered logarithmic potential. The outcome is the full spectrum of eigenfrequencies and their 
conjugate normal modes for prescribed azimuthal wavenumbers. The results show that thc funda- 
mental bar mode is isolated in the frequency space while spiral modes belong to discrete tamilies that 
bifurcate from the continuous tamily of van Kampen modes. The population of spiral modes in the 
bifurcating family increases by cooling the disk and declines by increasing the fraction of dark to 
luminous matter. It is shown that the variety of unstable modes is controlled by the shape of the dark 
matter density prohlc. 

Subject headings: stellar dynamics, instabilities, methods: analytical, galaxies: kinematics and dy- 
namics, galaxies: spiral, galaxies: structure 



1. INTRODUCTION 

Dynamics of self-gravitating stellar systems and 
plasma Auids are governed by the collisionles s Boltzmann 
equation (CBE) i|i3innev fc Tremaine1ll987[) . Einding a 
general solution of the CBE has been a challenging prob- 
lem in various disciplincs of physical sciences. Due to 
existing dithculties of the general problem, hnding a so- 
lution to the linearized CBE becam e the center of at trac- 
tion in the twentieth century when iLandau I (|1946l ) and 
van Kampen (1955) discovered the normal modes of coUi- 
sionless ensembles. Later in 1970's, Kalnais (1971,1977) 
developed a matrix theory that was capable of computing 
normal modes of stellar systems through solving a non- 
linear eigenvalue problem. His theory remained as the 
only analytical perturbation theory used by the commu- 
n ity of gala ct ic dyn amicists over the past three decades. 

iKalnaisl (|1977f ) assumed an exponential form 
exp(— iwt) for the time-dependent part of physical 
quantities where i = \/— 1, and solved the linearized 
CBE for the perturbed distribution function (DF) /i in 
terms of the perturbed potential Vi. After expanding 
the potential and density tunctions in terms of bi- 
orthogonal basis sets in the conhguration space, he used 
the weighted residual form of the fundamental equation 



fidvdx = Yiidx, 



(1) 



to obtain a nonlinear eigenvalue problem for thc complex 
eigenfrequency uj. For self-consistent perturbations thc 
surtace density Si = J fidv is related to Vi through 
Poisson's integral, and the symbols dv and dx denote 
th e elern e nts o f velocity and position vectors. 

IZang I (|l976f ) used Kalna.is' s theory to c ompute the 
modes of the isothermal disk (|Mestel II1963I ). which has 
the astrophysically important property of a flat rotation 
curye. His analysis was then extended bv lEvans fc Readl 
()1998al lbf) to general scale-free disks with arbitrary cusp 
slopes. Application of Kalnajs's theory to soft-centered 
models of stellar disks has been mai nly focused on 
the isochrone and Kuzmin-Toomre disks (|Kalnais Ill978t 



iHunter |[l99llPichon fc Cannon Ill997f) . A disk with ex- 
ponential light prohle and an approximately flat rotation 
curve was also investigated by Yauterin & Dcjonghe] 
(1996). More recently Jalali fc Hunter (2005a, herealter 
JH) gave new results for soft-centered models of stellar 
disks. They showed the importance of a boundary inte- 
gral in the modal properties of unidirectional disks and 
computed a tundamental bar mode and a secondary spi- 
ral mode for the isochrone, Kuzmin-Toomre and a newly 
introduced family of cored exponential disks. JH also ex- 
tended Kalnajs's first order perturbation theory to the 
second order, and illustrated energy and angular momen- 
tum content of different Eourier components. Their bar 
charts showed that only a few number of expansion terms 
in the radial angle govern the perturbed dynamics. 

Implementation of Kalnajs's (1977) theory, however, 
has some technical problems due to the nonlinear depen- 
dency of his matrix equations on to. Most computational 
methods that deal with nonlinear eigenvalue problems 
are iterative. They start with an initial guess of the 
solution and continue with a search scheme in the fre- 
quency space. Newton's method is perhaps the most ef- 
Acient technique that guarantees a quadratic convergence 
should the initial guess be close enough to the solution. 
The key issue in success of any iterative scheme is the 
attracting or repelling nature of an eigenvalue. It is ob- 
vious that only attracting eigenvalues can be captured 
by iterative methods while we have no priori knowledge 
of their basins of attraction in order to make our initial 
guess. The mentioned computational difflculties make 
it a formidable task to explore all normal modes of a 
stellar system, which include growing modes as well as 
stationary van Kampen modes. Moreover, it is not easy 
to develop a general nonlinear theory based on Kalnajs's 
method for studying the inte raction of modes. 

lPolvachenko I (|2004 I2005D introduced an alternative 
method for the normal mode calculation of stellar disks 
whose outcome was a linear eigenvalue problem for u>. 
His method is capable of hnding all eigenmodes of a stel- 
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lar disk should one use fine grids in the action space. 
Polyachenko's method is somehow costly because it re- 
sults in a large hnear system of equations to assure point- 
wise convergence in the action space and mean conver- 
gence of Fourier expansions in the space of the radial 
angle. Extension of his method to nonUnear regime 
is another challenging problem yet to be investigated. 
Tremaine ' ('2005') has also followed an approach similar to 
Polvachcnko (2005) and studied the instability of stellar 
disks surrounding massive objects. His eigenvalue equa- 
tions involve action variables, and practically, need to be 
solved over a discretized grid in the space of actions. 

Recent deyelop n ients in Auid me c hanics 
(iDoering fc Gibbon I [19951: lMattinglv fc Sinai I [l999h 
inspired me to formulate the dynamics of stellar systems 
in a new framework, which is capable of solving the CBE 
not only in the linear regime, but also in its full nonlinear 
form. The method systematically searches for smooth 
solutions of the CBE by expanding the perturbed DF 
using Eourier series of angle variables and an appropriate 
set of trial tunctions in the space of actions. CoelRcients 
of expansion are unknown time-dependent amplitude 
tunctions whose evolution egu ations are obtai ned by 
the Petrov-Galerkin projection ()Finlavson Ill972f) of the 
CBE. That is indeed taking the weighted residual form 
of the CBE by integration over the action-angle space 
and deriving a system of nonlinear ordinary differential 
equations (ODEs) for the amplitude functions. The 
associated hrst order system of ODEs leads to a linear 
eigenvalue problem, which is solved using standard 
numerical methods. 

In this paper I present my new method and apply it 
to explore the modal properties of a model galaxy. In 
a second paper, I will address the nonlinear evolution of 
modes and wave interactions. The paper is organized as 
toUows. In sections [2] and [3l I use the Petrov-Galcrkin 
method to project the CBE to a system of ODEs in the 
time domain and derive a system of linear eigenvalue 
equations. Section |4| presents the eigenfrequency spectra 
and their corresponding mode shapes of the cored expo- 
nential disk of JH. Thc stars of this model move in the 
field of a soft-centered logarithmic potential. I also inves- 
tigate the effect of physical parameters of the equilibrium 
model on the modal content. In section[5l I discuss on 
the nature of a self-gravitating mode and compare the 
pertormance of my method with other theories. Some 
tundamental achievements of this work are summarized 
in section[6l 

2. NONLINEAR THEORY 

I use the usual polar coordinates x — {R, (f)) and as- 
sume that the temporal evolution of the DF and gravita- 
tional potential starts from an axisymmetric equilibrium 
state described by fo{x,p) and Vo{R) so that 

f{x,p,t)=fo{x,p) + fi{x,p,t), (2) 
V{x,p,t) = Vo{R) + Vi{x,p,t). (3) 

Here p = {pr^p^,) is the momentum vector conjugate to 
X = (i?, (/)). Motion of stars in the equilibrium state is 
governed by the zeroth order Hamiltonian 

'Ho = UpR + ^]+MR)- (4) 



For bounded orbits, R and (f) become librating and ro- 
tating, respectively. One can therefore describe the dy- 
namics using the action variables J = {Jr,, J^), 

Jr= PRdR, J^ = (j) P4,d(f) = p^, (5) 

and their conjugate angles O = {6r,9^). A transtor- 
mation {x,p) (O, J) leaves the Hamiltonian Tio as a 
tunction of actions only, 7io(«/), and therctorc, the phase 
space Sows of the equilibrium state lie on a two dimen- 
sional torus J = c with c being a constant 2-vector. An 
action-angle transtormation can locally be found for any 
bounded regular orbit, but it is a global transformation 
if only one orbit tamily occupies the phase space. The 
axisymmetric potential Vo{R) supports only rosette or- 
bits. Radial and circular orbits are the limiting cases of 
rosette orbits with = and Jr = 0, respectively. By 
representing / in terms of the action-angle variables, the 
CBE reads 

^ + [./,W]=0, (6) 

where [, ] denotes a Poisson bracket taken ov er the action- 
anglc spacc. According to Jeans theorem peans I [19151 : 
[Lyndcn-BcU 1962) /o depends on the phase space co- 
ordinates through the integrals of motion, which are 
the actions in the present tormulation, and one obtains 
[/o,Wo] = 0. Subsequently, equation ^ may be rewrit- 
ten as 

^ = -[/i,Ho]-[/o,Hi]-[/i,7^i], (7) 

where Tii is the perturbed Hamiltonian. A dark matter 
halo contributes both to TLo and to TLi if it is live, i.e., if it 
exchanges momentum/energy with the luminous stellar 
component. In this paper I conhne myselt to a rigid halo 
that only contributes to TLo through Vo and assume that 
TLi = Vi is the perturbed potential due to self-gravity. 

Let me expand /i and Vi in Eourier series of angle 
variables and write 

oo oo 

fi{Q,J,t)= ^df'(t)$™'(J)e'("^*+'««\ (8) 

m.l— — oo j—0 

OG OG 

Vi{e,J,t)= ^6™'(t)*;"'(J)e'('"«*+'^«), (9) 

m,/— — oo j—0 

where $™'(J) and ^'"'(J) are some trial tunctions in 
the space of action variables, and c?™'(t) and fo™'(i) are 
timc-dependent amplitude functions. On thc other hand, 
one can expand Vi and its corresponding surtace density 
Ei in the conhguration space as 

oc oc 

^,{R,c^,t)=Y E^w^!'"'^^)^'"^ (10) 

m— — QO j—0 

oc oc 

v{R,c^,t)=Y E<wv^l™'(^)^'"'- (11) 

m— — oc j—0 

Here Vi"'(^) and ^'"'(i?) are bi-orthogonal potential- 
surtace density pairs that satisty the relation 

oo 

2n j i}^''\R)aY\R)RdR = Dj{m)5ra,m'5,,y, (12) 
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where 6m,m' is the Kronecker delta and Dj{m) are some 
constants. It is remarked that the real parts of Ei, Vi 
an d /i describe physical s olutions. In this work I utihze 
the lClutton-Brock l ()1972[ ) lunctions 



|m|. 



(13) 



|m| 



2|mH-2j + l\ /1-e 



3/2 



that yield (|Aoki fc Ive II1978I : iHunterlllOSl 

(2|m|+j)! 



P;"I(C), (14) 



Dj(m) = Dj(—m) 



2bjl 



(15) 



(0 s-^s associated Legendre lunctions with i = |m| + 
j. Clutton-Brock functions have a length scale b, which 
makes them suitable for reproducing the potential and 
surlace density of soft-centered models. The choice of 
this parameter is an important step in the calculation of 
normal modes. I will discuss on this issue later in 21 

Equating © and pTjl . multiplying both sides of the 
resulting equation by eTip[—i(l0ii + m9^)] and integrating 
over the {Oji, 9^)-space, lead to (see also Kalnajs 1977 and 
Tremaine & Weinberg 1984) 

oc oo 

j2 6f (i)*f ( j) =E (-^)' (16) 

j=0 j=0 



1 



*J"(J) = ^ / 4"'{R)cos[l9R + m{9^-cb)]d9RXl7) 



where are thc Eourier coethcients of the basis poten- 
tial tunctions in the conhguration space. The trial func- 
tions used in the expansion of Vi in the action-angle space 
are not necessarily identical to VE'™'. However, subse- 
quent mathematical derivations are greatly simphhed by 
setting ^f'{J) = Wf{J), which imphes bf\t) = af{t). 
To build a relation between (i™'(t) and a™{t), I use the 
fundamental equation 



fi{e,J,t)dJde = Si(i?, 0, t)RdRd(j), 



(18) 



where dJde is the volume of an inhnitesimal phase space 
element. On substituting ([5]) and (fTUl) in (fTS)) . multiply- 

ing both sides of the resulting equation by V'j"'' (^)e^'™'^ 
and integrating, one obtains 



<{t)- 



47r^ 



D,{m) 



EE 

i= — oo p=0 



A ml 

jP ~ 



J ^f{j)^;''{j)dj, 



(19) 



(20) 



which is inserted in ^ to represent Vi in terms of the 
amplitude tunctions d^^{t) as 



oc oo 



47r^ 



DAm) 

Tn,l,k— — oo j,p—0 



A™'=^'f'(J)d™'=(t)e'(™''*+'^«). 

(21) 

It would be computationally favorable to coUect d™' {t) 
in a single vector z{t) = {z„{t)} by dehning a map 



{m,l,j) — !■ n. In practice the inhnite sums in ([8]) 
are truncated and approximated by hnite sums so that 

-^max < l < 'max, -"lmax < m < m^a^ aud < j < 

imax- For 1 < n < nmax, a simplc map between indices 
will be 

n={m + mmax) (2Zmax + 1) (jmax + 1) 
+ {l+ ?max) {j max 

+ l)+j + l, (22) 

(2m,nax + 1) (2Zmax + 1) (jmax + 1)- (23) 

One can now use ([8]) and ((2T|) in ([7]) and apply the Petrov- 
Galerkin method to construct the weighted residual form 
of the CBE. That is to multiply ([7|) by some weighting 
tunctions W™'(0, J) and to integrate the identity over 
the action-angle space. The outcome is the following 
system of nonlinear ODEs 



dt ^ ■ 

q=\ 



pqZq ^ ^ Bpqj^ZqZ-p^ 



p = 1,2, - ■ ■ ,n„ 



for the amplitude functions z„(i) 



jml 



{t). 



(24) 

The elements 



of Apq and Bpgr have been determined in Appendix A. 



Each equation in ([241) is the projection of the CBE on 
a subspace spanncd by a weighting function. Therefore, 
the left hand side of ([24]) is the projection of dfi/dt, 
the summation over hrst order terms is the projection 
of — [/ij^o] — [/o,^i], and the second order terms are 
the projections of — [/i,Hi]. The second order terms 
of amplitude functions, characterized by Bp^r, show the 
interaction of modes in both the radial and azimuthal 
directions. 

Distribution of angular momentum between different 
Pourier components provides useful information of the 
disk dynamics. I compute the rate of change of the total 
angular momentum C using (see Appendix B in JH) 

d£ _ 1 

dt ~ "4. 



(/i + /i) ^ (^i + ^^i) dJde, (25) 



where a bar denotes complex conjugate. Substituting 
and (PT|) in (PS)) and evaluating thc integrals, yield 



dC 
'dt 



\.l=~oo j,p=0 



J2\m 4-™)(i)+^(t) A™'<(t) 



t^Tpdf{t) 



(26) 



Dehne af{t) = uf{t) + 'ivf{t) with u™(i) and y^p^t) 
being real lunctions of time. According to identity (|19p . 
one may lurther simplity equation (I26p to 



dL 

Itt 



m — — OQ 



(27) 



im(i)--f E^j("^) 
J=0 

X [uf{t)v-"\t)+uj"\t)vj^{t)]. (28) 

The share of the mth mode from dC/dt is thus deter- 
mined by Lm{t). As one could anticipate for an isolated 
stellar disk, dC/dt vanishcs and thc total angular mo- 
mentum remains constant because the terms Lm{t) and 
L~m{t) cancel each other in (|?71) . and Lo(0 is annuUed 
by the tactor m, in ((28| . 
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2.1. Trial and Weighting Punctions 

Choosing the trial functions $™'(J) is thc most dcl- 
icate step in the reduction of the CBE to a system of 
ODEs. One possible way is to set dfi/dt = in ([7]) and 
solve the first order equation 

[/i,7^o] + [/o,Hi] =0, 

for /i. Substituting © and ® in ([21]) gives 



df^f {J)^bf gf ( J)^™' ( J) , 



- dfo 



3 

dfa 



where 



dHo 



, ^4>{J) 



dHo 



Equation ((30| suggests to choose 

$"'(J) = g^^^J^-^^iJ) 



(29) 

(30) 
(31) 

(32) 
(33) 



as the trial lunctions (in the space of actions) for an 
unsteady /i(0, J,t). These functions have integrable 
singularities for resonant orbits with I^Ir + mO,^ — 0. 
For unidirectional disks with only prograde orbits, they 
also include a term with the Dirac delta function S{ J^) 
(see JH). One should therefore avoid the partial deriva- 
tives of with respect to the actions by evaluating the 
weighted residual form of [/i, Tii] through integration by 
parts (Appendix A). 

For deriving the relation between a™ and d™' in ((T9)) . 
the fundamental equation (|18p was multiplied by the 
complex conjugates of the basis functions used in the 
expansion of Vi {R, 4», t) . One may follow a similar ap- 
proach for obtaining the weighted residual form of the 
CBE and set 



VK™'(e, J) = ^f\j)e 



(34) 



which are the complex conjugates of the basis tunctions 
used in the expansion of Vi{Q, J ,t) in the action-angle 
space. The trial and weighting lunctions introduced as 
above, are not orthogonal but they result in a simple 
form for the linear part of the reduced CBE as I explain 



3. LINEAR THEORY 

In a first order perturbation analysis, the second order 
terms of the amplitude functions are ignored. The evo- 
lution of modes is then governed by the linear parts of 
(|ATo1) as 

iM ■ ■^^{t) = C ■ z{t). (35) 



dt 



A general solution of ([55]) has the form z{t) = 
which leads to the loUowing linear eigenvalue problem 



e Zq, 



C{m) ■ Zo = ujM{m) ■ Zc 



(36) 



for a prescribed azimuthal wavenumber m. Operating 
M"^ on (IMl) yields 



A{m) ■ Zq = ujZo, 



(37) 



where j4 is a general non-symmetric matrix. A reduc- 
tion to Hessenberg form followed by the QR algorithm 



(|Press et al. II2001D gives all real and complex eigenval- 
ues. Real eigenvalues correspond to van Kampen modes 
and complex eigenvalues, which occur in conjugate pairs, 
give growing/damping modes. I utilize the method of 
singular value decomposition for hnding the eigenvectors 
and perform the decomposition A — loI = U ■ S ■ 
where I is the identity matrix and the diagonal matrix S 
is composed of the singular values 5^^ (j = 1, 2, • • • , n^a^y^). 
The column of V that corresponds to the smallest Sj is 
the eigenvector associated with lo. 

Calculation of C and M involves evaluation of some 
dehnite integrals in the action space. There will be two 
types of such integrals (instead of three) if one uses the 
trial lunctions dehned in (|33|) . Let me introduce the aux- 
iliary integral 



= dJ l 



OJr 



^f{J)^f{J), 



(38) 



7nl 



and apply the trial lunctions $ 
The elements of M and C are thus computed from 



Mpq=dl,l'A'^y 



Cri 



SwT: 



ml 

n' 



E 

/c=0 



4^^ 



Dk{m) 



q-ml A 'ml' 



(40) 



Both AJ^' and consist of boundary integrals when the 



'-jk 



unperturbed stellar disk is unidirectional with the DF 
/o( J) = H{J^)fQ {J). Here H is the Heaviside function. 
The boundary terms are 



A ml 

^^]k - 



dJh 



mfP{J)^f{J)^f{J) 
lilR{J) + mn^{J) 



,(41) 



dJR [mf^{J)^f{J)^t{J) 



J J*=o 



J*=o ■ 



(42) 



Dynamics of modes with different azimuthal wavenum- 
bers are decoupled in the linear regime and the ma- 
trix A is an odd function of the wavenumber m. i.e., 
A{—m) = —A{m). An immediate result of this prop- 
erty is aj"^{t) = aj-{t). Consequently, L„i{t) becomes 
equal to zero for all |m| > and each mode individually 
conserves the total angular momentum. 

The present theory has three major advantages over 
Kalnajs's lormulation. Firstly, all eigenmodes relevant to 
a prescribed azimuthal wavenumber are obtained at once 
with classical linear algebraic algorithms. This makes it 
possible to explore and classify all families of growing 
modes beside pure oscillatory van Kampen modes. Sec- 
ondly, the constituting integrals of the elements of Alpg , 
Cpq and Kpqr (Appendix [X| are regular at exact reso- 
nances when the condition lQfj + mQ^ — cj = holds. 
Einally, nonlinear interaction of modes, and the mass 
and angular momentum exchange between them, can be 
readily monitored by integrating the system of nonlin- 
ear ODEs given in ((M)) . In the proceeding section I will 
be concerned with the calculation and classihcation of 
modes in the linear regime. 

4. MODES OF THE CORED EXPONENTIAL DISK 

JH calculated barred and spiral modes of certain stellar 
disks for the wavenumber m = 2. Among the models 
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FlG. 1. — Eigenfrequency spectra of a corcd exponential disk witli vo = 1, Rc = 1 and {N,X,a) = (6, 1,0.42). Eigenfrequencies have 
been displayed for the azimuthal wavenumbers < m < 5. The results correspond to (max = 10 and jinax = 15 in the series expansion of 
the perturbed distribution function. 



studied in JH, the cored exponcntial disk with the surtacc 
dcnsity prohlc 



havc derivcd thc gravitational potcntial corresponding to 
I dcnote this potential by Vd{R)- Thc gradicnt 



Rc 



Ec(i?) = S,cxp(^-A^l + i?2/i?2,j , A = ^, (43) 

and embeddcd in thc held of thc soft-centered logarith- 
mic potential 



FH = ^[V„iR)-VDiR)] 



(45) 



VoiR)^ vl In ^Jl + RyRl, (44) 

is a viable model that resemblcs most teatures of realis- 
tic spirals. Here Rc is the core radius, Rd is the length 
scale of the exponcntial dccay, and is a dcnsity scal- 
ing tactor. The velocity of circular orbits in this niodcl 
rises from zero at the galactic center and approaches to 
thc constant value vq in outc r rcgions whcrc the light 
prohlc falls off cxponentially. iJalali fc Huntcr I (|2005bD 



will givc thc gravitational forcc of a spherical dark mat- 
tcr componcnt, computed inside the galactic disk. The 
density profilc of the dark component, pn, can thcn 
bc determincd using Fh- The positivencss of pn im- 
poses some restrictions on the physical valucs of A and 
a = GT,sRd/vq as Pigure 5 in JH shows. Por a given A, 
a cannot cxceed a critical valuc «cr- Thc parametcr A 
dctermines the shape of thc dark mattcr dcnsity profile. 
A model with A = 1 and a = a^r is maximal in the region 
whcrc the rotation curve is rising. i.c, there is no dark 
mattcr in that rcgion. Models with A > 1 and a = a^r 
are still maximal but only in the vicinity of the ccnter for 
R < Rd- In such models the rotational velocity of stars 
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B1 (b) 




0.2 0.4 0.6 0.8 



COr 



FlG. 2. — (a) Zoomed eigenfrequency spectra (filled squares) of a cored exponentiaI disk witli {N,X,a) = (6,1,0.42) for the azimutlial 
wavenumber m = 2. The (isolated) bar mode has been labeled Bl. The most prominent growing modes belong to a discrete family 
that bifurcates from a van Kampen mode. The members of this family, labeled as Sl, S2, S3,- -, have spiral patterns. Circles show 
the eigenfrequencies of the tundamental and secondary modes calculated using Kalnajs's (1977) method (see Table 4 in JH). (6) The 
eigenfrequency loci of the same model of panel a as a is increased continuously from 0.2 to 0.42. A fcw sample eigenfrequencies have been 
displayed on each locus. 

error threshold of 1% the program terminates when 
('max, jmax) = (10, 15), which gives a size of 336 x 336 for 
the matrix A. In such a circumstance, out of 336 eigen- 
frequencies of A (for each wavenumber m), less than 15 
pair have non-zero growth rates (w/ 0). Further in- 
creasing of /max and jmax does not aher the number and 
location of complex eigenfrequencies in the w-plane. This 
shows that unstable modes do not constitute a continu- 
ous family. 

Figure [1] displays the eigenfrequency spectra for the 
azimuthal wavenumbers < m < 5. Eigenfrequen- 
cies on the real axis are oscillatory van Kampen modes. 
Their calc ulation reguires eva luation of Cauchy's princi- 
pal value ()Vandervoort 1120031 ) if one uses Kalnajs's first 
order theory. In the present formalism, van Kampen 
modes are found together with growing modes without 
any special treatment. More van Kampen modes are ob- 
tainable by increasing the truncation limit /niax of Fourier 
terms in the ^/{-direction. Toomre's Q is marginally 
greater than 1 for the model (Figure 7& in JH), and there- 
fore, one could expect that the disk is stable for m = 
excitations (see top-left panel in Figure[T]). The model is 
highly unstable for m > excitations although the av- 
erage growth rate of unstable modes decreases for larger 
wavenumbers. It is evident that either unstable modes 
are isolated or they are groupcd in discrete Jamilies. De- 
pending on the wavenumber, there may be one or more 
discrete families. The most prominent family bifurcates 
from van Kampen modes. Members of this family have 
spiral patterns with multiple peaks in their (i?) func- 
tions. The (global) fastest growing mode belongs to the 
spectrum of m = 2. That is the bar mode of a two- 
member unstable family. 

The length scale of Clutton-Brock functions has been 
set to & = 1.5 for TO = 2 and 6 = 2 for other wavenumbers. 
Changing this length scale slightly displaces the eigen- 
frequencies although the spectrum maintains its global 
pattern. Large values of b lead to a better computa- 
tion accuracy of extensive modes (with smaller pattern 
speeds), while compact bar modes show a rapid conver- 
gence for small values of b. Moreover, the suitable value 
of b differs from one azimuthal wavenumber to another. 
Finding an optimum length scale that gives the best re- 
sults for all modes and wavenumbers is an open problcm 



due to dark matter {vh — ^/RFh) has a monotonically 
rising prohle. For A < 1, dark matter penetrates into the 
galactic center and its density profile becomes cuspy in 
the limit of A ^ 0. The role of the parameter a is to 
control the fraction of dark to luminous matter. Models 
with a <C «cr are dominated by dark matter. 

JH introduced a family of equilibrium DFs that re- 
produces S£)(i?) and depends on an integer constant N . 
This parameter controls the population of near-circular 
orbits an d the disk temperature: the parameter Q of 
iToomre 1 ([l964) decreases by increasing N . The DFs of 
JH have an isotropic part that determines the fraction 
of radial orbits. That isotropic part, which reconstructs 
the central density of the equilibrium state, shrinks to 
central regions of the galaxy as N increases. 

I apply my new method to the cored exponential disks 
of JH and calculate the spectrum of w = ujn -\- iw/. Sub- 
sequently, the eigenvector Zo is calculated from ([37|) and 
it is used in ^ to compute aj'(t) = e~'"*a™(0) and the 
perturbed density 

Ei(i?, (/., t) = e"^*P™(i?) cos [m<t)-ujRt+'d^{R)\ , (46) 

which is the real part of (fTO)). Pm{R) and 'dm{R) are the 
amplitude and phase functions of an m-fold circumfer- 
ential wave that travels with the angular velocity ojnjm. 
The factor e'^'* shows the exponential growth/decay of 
the wave amplitude. I normalize all length, velocity, and 
time variables to i?c and Rc/vo, respectively, and 
set G — Rc = vo = 1. 

I begin my case studies in ^4.11 with a near maximal 
disk of (A^, A, a)=(6, 1, 0.42) and compute its eigenfre- 
quency spectra for the wavenumbers < "t, < 5. I then 
classify unstable m = 2 modes of this model and invcsti- 
gate their evolution as the parameter a is varied. In ^ 34.21 
and §4.31 I study the behavior of unstable m = 2 waves 
as the parameters A and N are changed. The eigenfre- 
quency spectrum of a model with an inner cutout is also 
computed and discussed in M.4\ 



4.1. A Near Maximal Disk 

I pick up the hrst model from Table 4 of JH with 
(A^, A,q:) = (6,1,0.42) and start solving the eigensys- 
tem ((57)) with (/max,imax) = (2,4) and increase these 
limits until complex eigenfrequencies converge. For an 
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FlG. 3. — Mode shapes of the cored exponential disk for Rc = 1, t)o = 1 and {N,X,a) = (6,1,0.42). Panels have been labeled by the 
corresponding mode name. The contour plots show the positive part of Ei (i?, iji, 0). The contour levels range from 10% to 90% of the 
maximum of Si (i?, 0, 0) with increments of 10%. The panel below each mode shape shows the amplitude of wave patterns as dehncd in 
equation II46I I. 



yet to be investigated precisely. Por the cored exponen- 
tial disks with 0.5 < Rc,Rd < 2, working in the range 
1 < 5 < 2.5 gives reasonable results. 

For m — 2,1 have zoomed out and plotted in Pigure 
^Bp, the portion of the spectrum that contains growing 
modes. The first and second modes reported in Table 4 
of JH have been shown by circles in the same hgure. The 
most unstable mode (labeled as Bl) is a compact, rapidly 
rotating bar. The majority of unstable modes belong to 
a discrete spiral tamily that biturcates from a van Kam- 
pen mode with lo « 0.43. I have labeled these modes by 
Sl,- • - jSG. The number of density peaks along the spiral 
arms is proportional to the integer number in the mode 
name. Both B2 and G are double peaked spirals but I 
have classitied B2 as a bar mode, and coUected it with Bl 
in a two member tamily, for it takes a bar-hke structure 
when it is stabihzed by decreasing a. I classity mode 
G as an isolated mode because it does not behave simi- 



lar to either of S- or B-modes as the model parameters 
vary. There is another isolated mode in the spectrum, 
Cl, which exhibits a spiral pattern. By decreasing A, 
mode Cl joins a new tamily of spiral modes, which are 
accumulated near the galactic center (see §4.2|) . 

Reducing a inc reases the abun dance of dark matter 
and according to iToomre 1 (jl981f ) and JH the growth 
rate of modes should decrease. My calculations show 
that by reducing a, spiral modes are affected sooner and 
more effective than the bar mode, and they join to the 
stationary modes, one by one from the location of the 
bifurcation point until thc whole S-family disappears. 
This is a generic scenario for all A > 1 models regardless 
of the disk temperature controlled by N . Sohd hnes in 
Pigure [2j' show the eigenfrequency loci of a model with 
{N, A) — (6, 1) as a increases from 0.2 to 0.42. It is evi- 
dcnt that mode Bl is destabihzed through a pitchtork bi- 
turcation while the loci of modes B2 and Cl, and S-modes 
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FlG. 4. — (a) Eigenfrequency loci (solid lines) of cored exponential disks with {N,X) = (6,0.625) for 0.2 < a < 0.36. Large squares 
show the eigenfrequency spectrum of a model with a = 0.34, and circles show the eigenfrequencies of the fundamental and secondary 
modcs calculated using Kalnajs's method (see Table 4 in JH) for a = 0.34. (6) Same as panel a but for models with {N, A) = (6, 2) and 
0.25 < a < 0.6. Some sample eigenfrequencies have been demonstrated on each locus. 



exhibit a tangent bifurcation. AU modes except mode 
G are stable for a < 0.23. Surprisingly, mode G resists 
against stabilization even for very small values of a. This 
indicates that mode G is not characterized by the frac- 
tion of dark to luminous matter. In ! j4.2l and ! j4.3[ I will 
show that this mode is highly sensitive to the variations 
of A and N. According to my computations (e.g., Figure 
[5^)), by increasing a all S-modes are born at the same bi- 
turcation frequency ujs ~ 0.43, but mode B2 comes out 
from a van Kampen mode with lub ~ 0.83. This result 
completely rules out any skepticism that mode B2 is a 
member of S-family. 

Figure [3] displays the wave patterns and amplitude 
tunctions of modes Bl, B2, G, Sl, S2 and S3. It is seen 
that the patterns of S-modes rotate slower and become 
more extensive as the mode number increases. Mode S6, 
which is at the bifurcation point of the spiral tamily, has 
the largest extent. Eight wave packets of this mode are 
distributed by a phase shift of 90 degrees along major 
spiral arms. Mode G has at most two density peaks on 
its major spiral arms but the magnitude of its second 
peak increases as the disk is cooled. Modes Bl, B2, Sl, 
S2 and S3 are, respectively, analogous t o modes A B, C , 
E and F of a Gaussian disk explored in lToomrĕl (|1981[). 
There are three low-speed modes in Eigure 1 1 of iToomre I 
([l98l5 that have not been labeled, but they are analogs 
of modes S4, S5 and S6. Mode G and Toomre's mode D 
also have some similarities but are of different origins (see 
None of them can be stabilized only by increasing 
the fraction of dark matter. 

Eigurc ^jp, shows that the tundamental mode obtained 
by JH coincides with mode Bl. The wave pattern of 
mode Bl (displayed in Figure|3]) is identical to the mode 
shape computed using Kalnajs's theory and demon- 
strated in Eigure 8 of JH. JH lound a secondary mode 
which lies between modes B2 and G. That mode is also a 
double-peaked spiral and it is not easy to identity its truc 
nature unless we investigate its evolution as the niodel 
parameters vary. By comparing Eigure lOa of JH with 
EigureHJ), one can see that both mode B2 and the sec- 
ondary mode of JH are destabilized through a tangent 
bifurcation while mode G has a different nature. The bi- 
turcation frequency of mode B2 that I hnd {ujb ~ 0.83) 
matches very well with the frequency of the stabilized 
secondary mode of JH (see Eigure lOa in JH but note 



that their vertical axis indicates Qp — uji^/2). There- 
fore, I conclude that the secondary mode of JH is indeed 
mode B2 although it seems to be closer to mode G. The 
existing discrepancy is due to different length scale of 
Clutton-Brock tunctions that JH have used for finding 
the secondary mode. By adjusting b one can improve the 
location of B2. However, this is an unnecessary attempt 
given the fact that mode B2 has already been identihed, 
and the computation accuracy of other eigenfrequencies 
has an impressive level. 

4.2. Yariations of X 

The parameter A controls the density profile of the dark 
matter component, specihcally near the galactic center. 
The traction of dark to luminous matter has its mini- 
mum value in marginal models with a ~ acr- I choose 
a marginal A < 1 model with (iV, A,») = (6,0.625,0.34), 
which has also been investigated by JH. EigureHla shows 
the portion of the spectrum that contains complex eigen- 
frequencies of this model. The spectrum has been com- 
puted for h — 1.5. Although m — 1 bar and spiral modes 
survive in this model, thcir pattern speeds and growth 
rates drop considerably. Mode G has been wiped out of 
existence by dark matter penetration into the center, and 
four unstable modes (Cl, C2, C3 and C4) have emerged 
that constitute a new family of spiral modes. They pop- 
ulate the central regions of the disk in most of A < 1 
models. Again, the location of eigenfrequencies obtained 
by JH have been marked by circles. The agreement be- 
tween the results of JH and the present study is very 
good and the variance is less than 2%. 

The population of spiral modes is changed by varying 
a, and the eigenfrequencics of unstable modes are al- 
tered significantly. Solid lines in EigurejUa show the loci 
of growing modes as a increases from 0.2 to a^r ~ 0.36. 
Similar to the previous A = 1 model, S-modes and mode 
Bl are destabilized through tangent and pitchtork bifur- 
cations, respectively. All C-modes are born by a pitch- 
fork biturcation although some minor modes of the same 
nature come and go as a varies. The loci of modes Bl 
and Sl (in Eigure IUa) are in harmony with the results 
of Kalnajs's method plotted in Pigure 106 of JH. It is 
noted that the locus of mode Bl steeply joins the real 
axis, well betore stabilizing the S-modes. This is how 
slowly growing spirals may dominate a stellar disk. 
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FlG. 5. — Same as Figure|3]but for a model with (A, a) = (0.625,0.34). The fastest growing spiral mode, Sl, has also been found by JH. 
Modes Cl, C2 and C3 have emerged due to dark matter presence at the galactic center. They shrink to central regions and grow slower as 
their pattern speed increases. 

There are no new families of growing modes in A > 1 
models. Dark matter in these models induces a rising 
rotation curve on the disk stars (see Pigure 6 in JH) and 
the population of S-modes dechnes. The growth rate of 
mode B2 increases proportional to A, but that of mode 
G faUs off although mode G is stiU robust against the 
yariations of a. Mode G has its maximum growth rate 
in A = 1 models, which suggests that it must be a self- 
gravitating response of the luminous matter that involves 
only the potential of the disk, Vd . Thc function {R) 
of mode Cl loses its minor peaks and becomes smoother 
as A increases. Figure|3^) shows the eigenfrequency loci 
of models with {N, X) = (6, 2) as a increases from 0.25 to 
acr ~ 0.6. The eigenfrequency loci of modes Cl and B2, 
and S-modes (as a varies) are similar to A = 1 models, 
but the locus of mode Bl loses its steepness and stretches 
towards small pattern speeds in an approximately linear 
form until it joins thc real w/{-axis. The biturcation fre- 
quency of modc B2 differs from S-modes and it is larger. 



Eigure [5] shows the wave patterns of modes Bl, B2, 
Sl, Cl, C2 and C3 for the model with (iV,A,a) = 
(6,0.625,0.34). The (isolated) mode Bl is stiU a single- 
peaked bar although its edge is more extensive as the 
Uat part of its P^iR) plot indicates. Mode Sl is a triple- 
peaked spiral (as before) and mode B2 is being stabilized 
(wb2 = 0.775 -I- 0.007i). It is seen that mode B2 has a 
bar-like structure, which justihes its classihcation as the 
secondary bar mode. The pattern of Sl and its P^iR) 
plot can be compared with Pigure 9 in JH. The agree- 
ment is quite satislactory. There is a remarkable differ- 
ence between the patterns of C- and S-modes although 
both families have spiral structures. In contrast to S- 
modes that become more extensive as their growth rate 
decays, C-modes are shrunk to central regions because 
their pattern speed increases. C-modes are also a bifur- 
cating family, but their biturcation point lies on large 
pattern speeds associated with the azimuthal frequency 
(Sl^) of central stars. 
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The parameters a and A are essentially controUing the 
fraction and density prohle of dark matter component, 
respectively. However, the radial velocity dispersion 
of the equilibrium state is also playing an important role 
in the perturbed dynamics. cth is an indicat or of thc ini- 
tial temperature of the disk. lEvans fc Read I ([l998b ) had 
already pointed out that the pitch angle of spiral patterns 
decreases when an ^ (see their Pigure 7). Apart from 
this morphological implication, can the variation in the 
disk temperature affect the modal content? In the follow- 
ing subsection, I trace the evolution of growing waves by 
changing the disk temperature and show that the popu- 
lation of S-modes is larger in rotationally supported, cold 
disks. 

4.3. Yariations of the Disk Temperature 

The parameter N of the DFs of JH controls the disk 
temperature by adjusting the size of thc isotropic core 
and the population of near circular orbits. As N in- 
creases, the streaming velocity (u^) approaches the ro- 
tational velocity of circular orbits and the stellar disk is 
cooled. Figure [6] displays the eigenfrequencies of previ- 
ous iX,a) = (1,0.42) and (A, a) = (0.625,0.34) models 
for 7V = 4 and iV = 8. The spectra for the intermediate 
value of = 6 have already been shown in Figures ^jp, 
and Hli. 

Increasing N gives birth to more S-modes while the 
bifurcation point of the family is preserved. As a new 
member is born at the bifurcation point, other members 
including mode Sl, are pushed away from the real axis on 
a curved path. This behavior is observed in both models 
but the branch of S-family in the model with (A, a) = 
(0.625, 0.34) stays closer to the real axis than the other 
model. The growth rates of C-modes increase remarkably 
as the disk is cooled. Despite mode Bl which rotates 
and grows laster in cold disks, mode B2 grows faster in 
warmer disks. Yariation in the disk temperature changes 
the eigenfrequency of modc G more effective than what 
a could, but nothing is more inAuential than the role of 
A. 

Another consequence of cooling the stellar disk is that 
m = waves are no longer stable. The parameter Q of 
Toomre was marginally greater than unity for A^ = 6 
models. For iV = 8, I find Q < l over an annular 
region because the plot of Q versus R exhibits a min- 
imum at some hnite radius (e.g., Figure 7 of JH). For 
instance, I find three growing m = modes for the 
model {N,X,a) = (8,1,0.42). They correspond to pure 
complex eigenfrequencies uji = 0.621i, uj^ = 0.494i and 
W3 = 0.238i. Mode shapes have (obviously) ringed struc- 
tures but the number of rings, which is identical to the 
number of peaks of Pq{R), depends on the growth rate. 
The modes associated with uji, loi and loz have three, 
four and five rings, respectively. Ring modes are very 
sensitive to the variations of model parameters and they 
are suppressed by decreasing A and a. 

4.4. The Ejject of an Inner Cutout 

In order to simulate an immobile bulge, which does not 
respond to density perturbations, JH utilized an inner 
cutout function of thc form 

-ffcut = l-e-(^*/^°)', (47) 
where Lq is an angular momentum scale. Multiply- 
ing Hcut by the self-consistent DF of the equilibrium 



state, prohibits the stars with < Lq from participat- 
ing in the perturbed dynamics. Consequently, incoming 
waves are reAected at some hnite radius and the inner- 
most wave packets of multiple-peaked modes are dimin- 
ished. My calculations show that all S-modes survive in 
cutout models, mode G disappears, and the growth rate 
of mode B2 increases. The pattern speed of mode Bl 
is boosted so that the corotation resonance is destroyed, 
but its growth rate drops drastically. Figure [7] shows the 
eigenfrequency spectrum of a model with Lo = 0.1 and 
{N, X, a) = (6, 1, 0.42). Circles show the eigenfrequencies 
found by JH. Again, the agreement between the results 
of JH and the present work is very good. The reason 
that I have identified mode B2 as the second member of 
B-family, and not the most unstable S-mode, is that its 
P2{R) function has an evolved double-peaked structure 
(see Figure 11 in JH) and its locus versus a does not 
emerge from the same bifurcation frequency of S-modes. 
Disappearance of mode G in cutout models conhrms my 
earlier note that it is a self-gravitating mode. 

5. DISCUSSIONS 

There are similarities between mode G of this study 
and Toomre's (1981) mode D. Both of these modes resist 
against stabilization by increasing the fraction of dark to 
luminous matter and they have at most double peaks on 
their spiral arms. Nonetheless, these modes are not the 
same because mode G is amplihed through a feedback 
from the galactic center but Toomre's mode D has been 
identified as an edge mode. A question remains to be an- 
swered: why Toomre (1981) did not detect mode G and 
I do not find an edge mode? The most convincing expla- 
nation is that to excite a self-gravitating wave inside the 
core of the stellar component, the governing potential in 
that region should mainly come from the self-gravity of 
stars. This requirement is fulfilled in my A > 1 mod- 
els for R < Rd. However, the completely hat rotation 
curve imposed by Toomre (1981) nowhere follows the ro- 
tational velocity induced by the self-gravity of stars and 
it prohibits the Gaussian disk from developing a G-like 
mode. On the other hand, I don't find an edge mode 
because the density protile of the cored exponential disk 
does not decay as steep as the Gaussian disk to create 
an outer boundary at some hnite radius for reAecting the 
outgoing waves. 

Similar to the hrst order analysis of ^ Polyachcnko's 
(2005) approach results in the full spectrum of eigen- 
frequencies for a given azimuthal wavenumber. There 
are some difFerences between his method and the present 
formulation. Polyachenko directly uses Poisson's integral 
to establish a point-wise relation in the action space be- 
tween the Pourier componcnts of the perturbed DF and 
its self-consistent potential. Combination of equations 
(5) and (9) in his paper is analogous to equation ([2T|) in 
this paper. The main departure of the two theories is in 
the way that the linearized CBE is treated. Polyachenko 
forces a point-wise fulfillment of the CBE in the action 
space while the present method works with a weighted 
residual form of the CBE. 

A point-wise formulation poses a challenge for the nu- 
merical calculation of the eigenvalues and thcir conju- 
gate eigenvectors. According to the bar charts of JH, at 
least ten Eourier components (—3 < / < 6) are needed 
in the ^^j-direction to assure a credible convergence of 
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FlG. 6. — The evolution of the eigenfrequency spectra as the disk temperature drops from left (7V : 
bottom panels correspond to (A,o) = (1,0.42) and (A,a) = (0.625,0.34), respectively. 
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FlG. 7. — Eigenfrequency spectrum of a cutout model with Lo = 
0.1 and {N,X,a) = (6,1,0.42). Circles show the eigenfrequencies 
found by JH using Kalnajs's theory. 

/i in a typical soft-centered galaxy model. Therefore, 
if onc chooses a grid of x in the action space, 
Polyachenko's eigenvector F will have a dimension of 
10 X ria X Tio. Therelore, for a very coarse grid with 
na — 21 that Polyachenko uses, thc unknown eigenvec- 
tor will have a dimension of 4410. This number must 
be compared with the dimension of Zo in equation ()37|) . 
That is indeed nniax = 336 for the most accurate calcula- 
tions carried out by setting (^inax, jmax) = (10, 15) which 
means that 21 Eourier components in the 0/j-direction 
and 16 expansion terms in the i?-direction have been 
taken into account. Noting that the detinite integrals 
Ij^' and AJl! are independently evaluated over the ac- 
tion space with any desired accuracy, the present theory 
proves to be more eSicient for eigenmode calculation (in 



the linear regime) than other existing alternatives. 

The agreement between the results of this work and 
those of JH, who have used Kalnajs's method, is im- 
pressive. There is only a discrepancy in the results for 
a double-peaked spiral mode of A = 1 models. In fact, 
these models have two double-peaked modes, modes B2 
and G, and JH hnd mode B2. The origin of discrepancies 
was attributcd to the lcngth scale of Clutton-Brock func- 
tions, 6, which is a fixed number for the whole spectrum 
of a given azimuthal wavenumber. Provided that JH op- 
timized b for each growing mode that they calculated (see 
also i )4.ip . some minor deviations from the results of this 
paper are reasonable. In most cases the algorithm used 
by JH converges to mode Bl and the fastest rotating S- 
mode. They capture mode B2 only if its growth rate is 
large enough. Other modes remain unexplored because 
Newton's method needs an initial guess of lu, which has 
a little chance to be in the basin of attraction of the 
other members of S-family. The separation of eigentre- 
quencies near the bifurcation point of S-modes is very 
small and one could anticipate complex boundaries for 
the basins of attraction of these eigenfrequencies. Thus, 
there is no guarantee that successive Newton's iterations 
keep an estimated eigenfrequency on the same basin that 
it was initially. Nevertheless, in Kalnajs's formulation, a 
systematic search for all growing modes i s possible by 
introducing the math ematical eigenvalue ()Zang I Il976t 
lEvans fc Read I fT998bf l and investigating its loci as the 
pattern speed and growth rate vary. 

6. CONCLUSIONS 

After three decades of Kalnajs's (1977) publication, 
it was not known exactly whether growing modes of 
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stellar systems appear as distinct roots in the eigenfre- 
quency space or they belong to continuous famiHes as 
van Kampen modes do. In this paper, I attempted to 
answer this question using the Galerkin projection of the 
CBE and unveiled the full eigenfrequency spectrum of 
a stellar disk. I showed that similar to gaseous disks 
(|Asghari fc Jalali"! I2006D . majority of growing modes 
emerge as discrete /amilies through a bifurcation from 
stationary modes. There are some exceptions for this 
rule, the most important of which are the isolated bar 
and G modes. 

The model that I used to test my method allows for 
dark matter presence as a spherical component, whose 
potcntial inside the galactic disk contributes to the ro- 
tational velocity of stars. By varying the parameters 
of the model, and investigating the eigenfrequencies and 
their associated mode shapes, I showed that it is not the 
traction of dark to luminous matter that controls the va- 
riety of growing modes. What determines that variety is 
indeed the shape of the dark mattcr density prohle con- 
troUed by the parameter A = Rc/Rd- My survey in the 
parameter space revealed that the concentration of dark 
matter in the galactic center (A < 1) destroys mode G 
and weakens the growth of B-modes substantially. Emer- 
gence of spiral G-modes that accumulate near the galac- 



tic center is another remarkable consequence of dark mat- 
ter presence in central regions of a cored stellar disk. 

Although the solution of the Galerkin system showed 
a credible convergence of the series expansions, the ex- 
istence of strong solutions for thc GBE, in its full non- 
linear form, is still an open problem. It has been known 
for years th at van Kampen modes make a complete set 
(jCase Ill959( ). and therefore, they may be used for a series 
representation of stationary oscillations. But there is not 
a mathematical proof for the completeness of the discrete 
tamilies of growing modes. In other words, whether an 
observed galaxy can be assembled using the modes of a 
linear eigensystem, requires turther analysis. 

In the second part of this study, I will investigate the 
mechanisms of wave interactions in the nonlinear regime 
and will probe the mass and angular momentum transter 
between waves of different Eourier numbers. 
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APPENDIX 



WEIGHTED RESIDUAL EORM OF THE COLLISIONLESS BOLTZMANN EQUATION 

Let me define a nonlinear operator A and denote m'^^-' as the £th prolongation (|01ver Ill993f) of the physical quantity 
u in the domain of independent variables. Assume a (nonlinear) partial differential equation 



A^yl-^^^^t^ = 0, 



(Al) 



and its associated initial and boundary conditions tha t govern the evol ution of u{x,t) in the domain of the spatial 
variable x and the time t. A weighted residual method (jFinlavson Ill972f l attempts to find an approximate solution of 
the form 

u{x,t) = ^ a},{t)(p},{x), (A2) 

k=\ 

through dctcrmining thc time-dependent tunctions ak(t) for a given set of trial (basis) functions ipk{x). The trial 
tunctions should satisty the boundary conditions and be linearly independent. Using (|A2[1 and taking the inner 
product of (lAip by some weighting functions Wk'{x), yield the determining equations of ak{t) as 

{A,Wk') = jAWk'dx = 0, fc' = l,2,---,A:,„ax. (A3) 

There are several procedures for choosing Wk' {x) and each procedure has its own name. The method with Wk' = (pk' is 
called the Bubnov-Galerkin, or simply the Galerkin method. The Pctrov-Galcrkin mcthod is associated with Wk' ^ (pk' ■ 
The well-known collocation method uses Dirac's delta tunctions for thc wcighting purposc. Thcrc is an alternative 
interpretation for the inner product {A, Wy) = 0. That is projecting the equation ^ = on a subspace spanned by 
the weighting tunction Wy ■ Theretore, equation (IA3|) is often called the Galerkin projection of (jAip . In what loUows, 
I use the Petrov-Galerkin method and construct thc wcighted residual form of the GBE. 
Assume the tunctions U{Q, J) and V{Q, J), and dehne their inner product over the action-angle space as 



{u,v) = j J u{e,j)v{e,j)djdQ. 



Taking thc inncr product of the perturbed GBE by the wcighting tunctions Wp^O, J) = W^^{J)e 

= " {[h.no],wf) - {[h,ni],wf) - {[h,n\],wf) . 



(A4) 
gives 
(A5) 



Note that the GBE is the governing equation of the pcrturbed DF whose trial functions are <I>™'(J). With my choice 
of the weighting function (as above) I am following the Petrov-Galerkin method. On substituting ([8]) and ((2T|) in (jA5[) 
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and after some rearrangements of summations, one obtains 
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m' d' j' 



m' ,1' j 



'EEw<''(^)E 

m' ,1' j' k 



A™/' I dJU 



m' ,1' j' m" ^l" j 



Dk{m') 

m" ,{m—m')^j' ''{t)df'"{t)j2 



<Yfm'l ( T\\Ttm'l 



dJB. dJi 

47r2 



^'"'■'(J)*^ (J) 



Dk{m") 

/.j*f(-')*r''(-')(''j^+'"'^)*f"-'''y) 



\m L 



dJ^f{J)^f^'-'\j)l{l^l')^+m"^j^^-^'{J) 



i'l' , 



(A6) 



where -TOmax < m,m',m" < mraa.^, -Zmax < 1,1' .,1" < /max and < j,j',j",k < jmax- Using equation and 
carrying out the index mappings (m, l, j) — + p, {m' , l',j') — + q and {m" , l" ,j") — ^ r one may introduce thc arrays 



M J dJ-^f{J)^f'{J), 

Si^i' J dj {i'nR + m'n^)^f{j)^f'{j) 



(A7) 



E 

fc=0 



" 47r2 " 






Kf 


jdj[ 


_Dk{m')_ 



(A8) 



-,ml ( T \jf^m'l' 



47r^ 



Dk{m") 



\ 7n"l" 

J-^kj" 



dJ^J\J)^f^ {J) {l'±-+m'^Uf'"{J) 



- J dJ^f{J)^f^"{J) (^l'' 



dJR 



Consequently, equation (|A6[) takes the following matrix form 



' E = E ^PI^^i*) + E KpqrZq{t)Zr{t), Zp{t) = df{t), p = 1, 2, • • • , « 

g=l q=l 

Let the matrix M"^ = [Mp,^] bc the inverse of M = [Mpg] and left-muhiply (jAlOp by M"^ to get 



i^^p(i) = £ Ap,z,(i) + 2 BpqrZq{t)zr{t), A = M-1 • C, Bpqr = £ Mp^ii^, 



(A9) 



(AIO) 



(All) 



9=1 



q.r—l 



Evaluation of the integrands in (jA9p will be considerably simplihed if one avoids the partial derivatives of $™'(J) 
through integrating (|A9|1 by parts. That gives 



Kpqr — Sm' ^(m-m")Sl' ^(l_l") 

fc=0 



J dJ<^>f'{J)^f'"{J) (^l" 



Dk{m") 

d_ _d 

dJii dJt 
d 



dJ-^f{J)^f'{J) ( l^+m^ ) ^f"'{J) 

d 



dJR dJcj, 



*7'(j) 



(A12) 
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When all stars move on prograde orbits, the equilibrium DF takes the form fo{J) = H{Jff,)fQ [J) where H is the 
Heaviside function. Upon using ([33|) . this contributes a term including Dirac's delta tunction 5{J^) to the trial 
functions. Thus, the toUowing boundary terms 



K" ^ S 



fc=0 



Dkim") 



m"l" 
kj" 



dJii 



dJn 



m'f^{J)^f{J)^f^'{J) 

i'nR{j) + m'n^{j) 

m'fP{J)Wf'"{J)Wp'''{J) 

i'nR{j) + m'n^{j) 



d 



d 



Jd,=0 



dJr, 



(A13) 



must be added to Kpqr when the equihbrium disk is unidirectional. The partial derivatives of W^^{J) needed for 
equations (IA12P and ()A13P are calculated by differcntiating equation (I17|l partiaUy with respect to an action: 



TjT 



dR dJu 



cos[Z0/ 



m^^P{R)^ {O^ - 0) &m[ieR + m( 



d 



der- 



R, 



(A14) 



Jalah & Hunter (2005b) encountered these partial derivatives in their second order perturbation theory devised for 
computing the energy of eigenmodes. I adopt their tcchnique for calculating the quantities dR/dJi, and d {e^ — 4>) / dJu. 
The variables i?, (^^ — i^), and pn are regarded as tunctions of {Jr, J^, 0_r) because the action-angle transformation 
{x,p) {J,Q) is defined in the phase space of an axisymmetric state. From = dR/dt — {dR/ deji)Q.}i one may 
write 



d 
di 

Similarly, one obtains 

d_ 

dt 



dR 



dJ, 



dJ, 



d'R deR 
deRdJ^ ~dr 



d f dR 



'dJ, \deR 



d 



dj^ \n 



vr 



dvR vr dilR 



dJu ^R dJn 



{O^ - 0) 



d_ 

di 



dvR 



dJu 



dfl4, _ 
dJ, 

2Jp. 
i?3 



2J^ dR 



1 

dR 
dJ, 



J(j, 

R^ 



1 



d^R 
dJ, ' 



dflR 
dJ, ■ 



(A15) 

(A16) 
(A17) 



The set of three equations (jAlsp through (|A17p can be integrated along an orbit, and they provide the additional 
values needed to evaluate the partial derivatives (|A14p . Initial values are dvR/dJ^ = 9(0^ — (^/dJ^ = Oat0/; = i = O 
where R = iJmin because vr = B^ — <^ = Q for all orbits. However the initial i?inin values change with the actions, and 
initial values for the derivatives of R with respect to the actions are 



dR 
dJ^ 



i?3 n 



R 



R=B.„ 



-^min^o'(-^min) ^1 



dR 



dJd, 



i?min (i?niin^0 



J^) 



R=R„ 



-^min^'(^min) Ja, 



(A18) 



They are obtained by differentiating the zcroth order energy equation. 
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